Starting with the Schaefer 400-parcel, 7 network atlas.
#set variables for `collect_data`
if(fslong){
fname <- 'sea_rsfc_fslong_schaefer400x7.RDS'
data_dir <- '/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/SEA_REST/derivatives/fmriprep-20.1.1/xcpengine-default/'
subs <- sprintf('sub-10%02d', c(1:16,18:31))
sess <- sprintf('%02d', 1:10)
subses <- expand.grid(sess, subs)
files <- paste0(data_dir, '/', subses[,2], subses[,1], '/fcon/schaefer400x7/', subses[,2], subses[,1], '_schaefer400x7.net')
file_id <- paste0(subses[,2], '_', subses[,1])
} else {
fname <- 'sea_rsfc_schaefer400x7.RDS'
data_dir <- '/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/SEA_BIDS/derivatives/1.4.1-final/xcpengine-default'
subs <- sprintf('sub-10%02d', c(1:16,18:31))
sess <- sprintf('ses-%02d', 1:10)
subses <- expand.grid(sess, subs)
files <- paste0(data_dir, '/', subses[,2], '/', subses[,1], '/fcon/schaefer400x7/', subses[,2], '_', subses[,1], '_schaefer400x7.net')
file_id <- paste0(subses[,2], '_', subses[,1])
}
#This should be easy to generalize
#option to set directory and search path
#option to rad in label file
library(data.table)
library(tidyr)
if(is.na(ncores <- as.numeric(Sys.getenv('SLURM_CPUS_ON_NODE')))){
ncores <- 10
message('No environment variable specifying number of cores. Setting to ', ncores, '.')
}
setDTthreads(ncores)
if(file.exists(fname)){
adt_labels <- readRDS(fname)
schaefer_400_7_net_ids <- readRDS('schaefer_ids.RDS')
} else {
message('Creating file list...')
files_list <- lapply(files, function(x) ifelse(file.exists(x), x, ''))
message('Reading rsfc files...')
if(ncores > 1){
message('Using ', ncores, ' cores to read files.')
cl <- parallel::makePSOCKcluster(max(c(ncores - 1), 1))
adt_list <- unlist(parallel::parLapply(cl = cl, split(files, 1:ncores), function(part){
lapply(part, function(f){
if(file.exists(f)){
data.table::fread(f, skip = 2, col.names = c('node1', 'node2', 'r'))
} else {
data.table::data.table()
}
})
}), recursive = FALSE)
try(parallel::stopCluster(cl))
} else {
adt_list <- lapply(files, function(f){
if(file.exists(f)){
data.table::fread(f, skip = 2, col.names = c('node1', 'node2', 'r'))
} else {
data.table::data.table()
}
})
}
names(adt_list) <- file_id
message('Combining data, assigning labels...')
adt <- rbindlist(adt_list, idcol = 'id')
#https://github.com/ThomasYeoLab/CBIG/tree/master/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal/Parcellations/MNI
schaefer_400_7_net_ids <- fread('Schaefer2018_400Parcels_7Networks_order.txt',
col.names = c('node1', 'network1', 'V1', 'V2', 'V3', 'V4'))
schaefer_400_7_net_ids[, c('V1', 'V2', 'V3', 'V4') := NULL]
schaefer_400_7_net_ids <- schaefer_400_7_net_ids %>%
tidyr::extract(network1, c('network1', 'roi1'), '7Networks_([RL]H_.*)_(\\d+)')
setkey(adt, node1)
setkey(schaefer_400_7_net_ids, node1)
adt_labels1 <- schaefer_400_7_net_ids[adt]
setnames(schaefer_400_7_net_ids, c('node2', 'network2', 'roi2'))
setkey(adt_labels1, node2)
setkey(schaefer_400_7_net_ids, node2)
adt_labels <- schaefer_400_7_net_ids[adt_labels1]
adt_labels[, c('id', 'sess') := tstrsplit(id, '_', fixed = TRUE)]
message('Saving RDS file to: ', fname)
saveRDS(adt_labels, fname)
message('Saving Schaefer ID data table to schaefer_ids.rds')
saveRDS(schaefer_400_7_net_ids, 'schaefer_ids.RDS')
}
if(FALSE){
adt_labels_old <- readRDS('sea_rsfc_schaefer400x7.RDS')
missing_files_csv <- data.table(file = files, does_not_exist = files_list == '')
View(missing_files_csv[does_not_exist == TRUE])
readr::write_csv(data.table(file = files, does_not_exist = files_list == ''),
'~/sea_rsfc-missing_fcon_output.csv')
a <- data.table(file = files, does_not_exist = files_list == '')
a[, fcon_missing := lapply(gsub('(.*/fcon)/.*', '\\1', file), function(x) !file.exists(x))]
a[, missmatch := does_not_exist != fcon_missing ]
a[(missmatch), file]
}
Retain just within-network connectivity. Also, Fisher-z transform the correlations for analyses. One small detail that is important here is that we keep network connectivity for same-hemisphere parcels only.
sea_schaefer400_7_withinnet <- copy(adt_labels)
sub_regex <- '[LR]H_[A-Za-z]+_*(.*)'
net_regex <- '[LR]H_([A-Za-z]+)_*(?:.*)'
hem_regex <- '([LR]H)_[A-Za-z]+_*(?:.*)'
sea_schaefer400_7_withinnet_nets <- distinct(sea_schaefer400_7_withinnet, network1)
sea_schaefer400_7_withinnet_nets[, sub1 := gsub(sub_regex, '\\1', network1)]
sea_schaefer400_7_withinnet_nets[, net1 := gsub(net_regex, '\\1', network1)]
sea_schaefer400_7_withinnet_nets[, hem1 := gsub(hem_regex, '\\1', network1)]
setkey(sea_schaefer400_7_withinnet, network1)
setkey(sea_schaefer400_7_withinnet_nets, network1)
sea_schaefer400_7_withinnet <- sea_schaefer400_7_withinnet_nets[sea_schaefer400_7_withinnet]
setnames(sea_schaefer400_7_withinnet_nets, c('network2', 'sub2', 'net2', 'hem2'))
setkey(sea_schaefer400_7_withinnet, network2)
setkey(sea_schaefer400_7_withinnet_nets, network2)
sea_schaefer400_7_withinnet <- sea_schaefer400_7_withinnet_nets[sea_schaefer400_7_withinnet]
sea_schaefer400_7_withinnet <- sea_schaefer400_7_withinnet[net1 == net2 & hem1 == hem2]
sea_schaefer400_7_withinnet[, c('network1', 'network2') := NULL]
sea_schaefer400_7_withinnet[, z := atanh(r)]
sea_schaefer400_7_withinnet[, H := dplyr::case_when(hem1 == 'RH' ~ -1,
hem1 == 'LH' ~ 1)]
if(!dim(distinct(sea_schaefer400_7_withinnet, net1, net2, hem1, hem2))[[1]] == 14){
stop('Incorrect number of networks')
}
Using the number nodes are in each network allows computation of the expected number of rsfc features, which allows a comparison to the observed number of features to find instances where certain connections are missing.
#Check to make sure we're getting the number of connectivity features per network we expect
setkey(sea_schaefer400_7_withinnet_nets, 'network2')
setkey(schaefer_400_7_net_ids, 'network2')
nodes_per_net <- sea_schaefer400_7_withinnet_nets[schaefer_400_7_net_ids][, .N, by = c('net2', 'hem2')]
nodes_per_net[, expected_Nfeatures := N * (N - 1) / 2]
connectivity_feature_check <- sea_schaefer400_7_withinnet[, .(Nfeatures = .N), by = c('net2', 'hem2', 'id', 'sess')][nodes_per_net[, !'N'], on = c('net2', 'hem2')]
setnames(connectivity_feature_check, c('net2', 'hem2'), c('nework', 'hemisphere'))
dt_options <- list(rownames = FALSE,
filter = 'top',
class = 'cell-border stripe',
extensions = 'Buttons',
options = list(dom = 'Bfrtip', buttons = c('csv')))
do.call(DT::datatable,
c(list(connectivity_feature_check[Nfeatures != expected_Nfeatures],
caption = 'Mismatch of observed and expected number of rsfc features'),
dt_options))
do.call(DT::datatable,
c(list(connectivity_feature_check[Nfeatures != expected_Nfeatures][, .N, by = c('id', 'sess')],
caption = 'Number of mismatches for subject sessions'),
dt_options))
Density of functional connectivity (pearson correlation) for each participant, for each subject, overlayed on the density for all sessions and participants.
#Generate group-level density
agg_mean <- mean(sea_schaefer400_7_withinnet$z)
agg_sd <- sd(sea_schaefer400_7_withinnet$z)
scalez <- (sea_schaefer400_7_withinnet$z - agg_mean) / agg_sd
density_agg <- density(scalez)
sea_schaefer400_densplot_agg <- data.table(x = density_agg$x,
x_on_z = density_agg$x * agg_sd + agg_mean,
y = density_agg$y)
#Generate per-participant-session densities
sea_schaefer400_MSD <- sea_schaefer400_7_withinnet[, list(mean = mean(z),
sd = sd(z)),
by = c('id', 'sess')]
sea_schaefer400_densplot <- sea_schaefer400_MSD[sea_schaefer400_7_withinnet,
on = c('id', 'sess')]
sea_schaefer400_densplot[, scalez := (z - mean)/sd]
sea_schaefer400_densplot <-
sea_schaefer400_MSD[sea_schaefer400_densplot[, list(x = density(scalez)$x,
y = density(scalez)$y),
by = c('id', 'sess')],
on = c('id', 'sess')][, x_on_z := sd*x + mean]
u_ids <- unique(sea_schaefer400_densplot$id)
for(an_id in u_ids){
cat(paste0('\n\n## ', an_id, '\n\n'))
hplot <- ggplot(sea_schaefer400_densplot_agg, aes(x = tanh(x_on_z), y = y)) +
geom_ribbon(ymin = 0, aes(ymax = y, fill = 'Group', color = 'Group'),
alpha = .5) +
geom_ribbon(data = sea_schaefer400_densplot[id == an_id, ],
ymin = 0, aes(ymax = y, fill = 'ID', color = 'ID'),
alpha = .5) +
scale_fill_manual(aesthetics = c('color','fill'), breaks = c('Group', 'ID'), labels = c('Group', 'Participant'), values = apal[c(2,4)], name = 'Data') +
facet_wrap(~sess, ncol = 2) +
labs(x = 'correlation', y = 'density') +
coord_cartesian(xlim = c(-1, 1), ylim = c(0, .5)) +
jftheme
print(hplot)
}
We are interested in stability and variability for each network. For each participant, we have multiple sessions, and within each session, we have multiple parcel-parcel pairs that provide information about the within-network connectivity. We can use the fact that we have multiple indicators of within-network connectivity to estimate the between-person variability, as well as the within-person (session-to-session) variability as distinct from measurement error (we assume that each parcel-parcel functional connectivity estimate in a network is an estimate of that network's connectivity)*.
* This assumption means that we essentially treat differences in the FC between one pair and another as measurement error.
We can compute an ICC that describes both within-person and between-person variability by using a 3 level model. First, I subset the data for a single network. I then build an intercept-only model, allowing the intercept to vary by participant-ID, and by session within each ID. The intercept is effectively the mean across all intrahemispheric parcel-pairs.
\[ \begin{align} z &= \beta_{0jk} + \epsilon_{ijk} \\ \beta_{0jk} &= \pi_{00k} + \nu_{0jk} \\ \pi_{00k} &= \gamma_{000} + \upsilon_{00k} \end{align} \]
\[ \begin{align} z = \gamma_{000} + \upsilon_{00k} + \nu_{0jk} + \epsilon_{ijk} \end{align} \] Where \(z\) is the measure of functional connectivity, and \(i\) indexes observations within each session, \(j\), for each participant \(k\). This means that we are able to estimate variance in per-person deviations (\(\upsilon_{00k}\)) from the network mean-connectivity (\(\gamma_{000}\)), per-session-deviations (\(\nu_{0jk}\)) from those per-person deviations, and the residual not explained by variability over persons or person-sessions (\(\epsilon_{ijk}\)).
I want to make sure that these models are correct, so I'll compare the variance portions, and total variance, between the two models.
library(lme4)
networks <- unique(sea_schaefer400_7_withinnet$net1)
this_net <- networks[[2]]
this_net_dt <- sea_schaefer400_7_withinnet[net1 == this_net]
fslong_insert <- ifelse(fslong, '_fslong', '')
model_dir <- 'models'
test_model_list <- list(test_2l = list(fn = file.path(model_dir, paste0('test_noh', fslong_insert, '_2l.RDS')),
form = z ~ 1 + (1 | id)),
test_3l = list(fn = file.path(model_dir, paste0('test_noh', fslong_insert, '_3l.RDS')),
form = z ~ 1 + (1 | id/sess)))
cl <- parallel::makePSOCKcluster(max(c(ncores - 1), 1))
test_model_fits <- parallel::parLapply(cl = cl, test_model_list, function(model_list, d){
library(lme4)
if(!file.exists(model_list[['fn']])){
f <- model_list[['form']]
fit <- lmer(f, data = d)
saveRDS(fit, model_list[['fn']])
} else {
fit <- readRDS(model_list[['fn']])
}
return(fit)
}, d = this_net_dt)
try(parallel::stopCluster(cl))
three_level <- test_model_fits$test_3l
summary(three_level)
## Linear mixed model fit by REML ['lmerMod']
## Formula: z ~ 1 + (1 | id/sess)
## Data: d
##
## REML criterion at convergence: 182326.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7489 -0.6822 -0.0861 0.5939 6.3316
##
## Random effects:
## Groups Name Variance Std.Dev.
## sess:id (Intercept) 0.0045427 0.0674
## id (Intercept) 0.0005854 0.0242
## Residual 0.0813511 0.2852
## Number of obs: 550247, groups: sess:id, 282; id, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.185095 0.005989 30.9
lll_varcorr <- VarCorr(three_level)
report_df <- data.frame(
stat = c('id_var', 'sess_var', 'resid_var', 'total_var'),
three_level = c(lll_varcorr$id.1['(Intercept)','(Intercept)'],
lll_varcorr$sess.id.1['(Intercept)','(Intercept)'],
sigma(three_level)^2,
sum(unlist(lapply(lll_varcorr, diag))) + sigma(three_level)^2))
report_df$three_level_perc <-
sprintf('%.1f%%', report_df$three_level/report_df$three_level[[4]]*100)
report_df$three_level_RE_perc <-
c(sprintf('%.1f%%', report_df$three_level[1:2]/sum(report_df$three_level[1:2])*100), '-', '-')
knitr::kable(report_df, digits = 5)
| stat | three_level | three_level_perc | three_level_RE_perc |
|---|---|---|---|
| id_var | 0.08135 | 94.1% | 48.5% |
| sess_var | 0.08648 | 100.0% | 51.5% |
| resid_var | 0.08135 | 94.1% | - |
| total_var | 0.08648 | 100.0% | - |
networks <- unique(sea_schaefer400_7_withinnet$net1)
netlist <- lapply(networks, function(this_net){
return(list(network = this_net, d = sea_schaefer400_7_withinnet[net1 == this_net]))
})
cl <- parallel::makePSOCKcluster(max(c(ncores - 1, 1)))
model_fits <- parallel::parLapply(cl = cl, netlist, function(anet, fslong){
this_net <- anet[['network']]
this_net_dt <- anet[['d']]
model_dir <- 'models'
fslong_name <- ''
if(fslong)
fslong_name <- 'fslong-'
model_fit_fn <- file.path(model_dir, paste0('lmer_fit-3l-', fslong_name, this_net, '.RDS'))
library(lme4)
if(!file.exists(model_fit_fn)){
fit <- lmer(z ~ 1 + (1 | id/sess), data = this_net_dt)
saveRDS(fit, model_fit_fn)
} else {
fit <- readRDS(model_fit_fn)
}
return(fit)
}, fslong = fslong)
try(parallel::stopCluster(cl))
proportion_variance_tables <- lapply(model_fits, function(model_fit){
vc <- VarCorr(model_fit)
id_v <- vc$id['(Intercept)','(Intercept)']
idsess_v <- vc$`sess:id`['(Intercept)','(Intercept)']
s_2 <- sigma(model_fit)^2
total_RE <- sum(unlist(lapply(vc, diag)))
other_RE <- total_RE - id_v - idsess_v
total <- total_RE + s_2
rez <- data.frame(source = rep(c('ID', 'ID/Session', 'Other RE', 'Residual'),2),
out_of = factor(c(rep('total', 4), rep('RE', 4)), levels = c('total', 'RE')),
percent = c(c(id_v, idsess_v, other_RE, s_2)*100 / total,
c(id_v, idsess_v, other_RE, NA)*100 / total_RE))
if(other_RE == 0){
rez <- rez[rez$source != 'Other RE',]
}
return(rez)
})
The plots on the left, below, show model-expected functional connectivity pearson correlations for each participant, for each session, overlayed on the raw data (one point per parcel, per session, per participant). A line for each participant's model-expected mean across all sessions is overlayed on this. The right plot shows proportions of variance accounted for by the random effect estimates as a proportion of both total variance, and total random effects (RE) variance. The total RE variance includes individual means (ID), individual means per session (ID/Session), and the difference in connectivity between right and left hemispheres (we treat this as a nuisance variable, essentially).
library(patchwork)
plot_percents <- function(percent_table){
ggplot(percent_table, aes(y = percent, x = out_of, fill = source)) +
geom_col(position = 'stack') +
scale_fill_manual(breaks = c('ID', 'ID/Session', 'Other RE', 'Residual'),
values = apal[c(4,3,2,1)]) +
scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100), labels = paste0(c(0, 20, 40, 60, 80, 100), '%')) +
labs(x = 'Variance (denominator)', y = '', fill = 'Variance source') +
jftheme
}
plot_predictions <- function(amodel, point_alpha){
adt <- as.data.table(amodel@frame)
adt[, wave := factor(as.numeric(gsub('ses-(\\d+)', '\\1', sess)))]
newdata <- unique(adt, by = c('id', 'sess', 'wave'))[, c('H', 'z') := list(0, NULL)]
newdata$z_prime <- predict(amodel, newdata = newdata)
newdata$z_prime_id <- predict(amodel, newdata = newdata, re.form = ~(1|id))
ggplot(newdata, aes(x = wave, y = tanh(z_prime), group = id)) +
#geom_violin(data = adt, aes(group = NULL, y = tanh(z)), size = 0, alpha = .5, fill = apal[[1]]) +
geom_point(data = adt, aes(group = NULL, y = tanh(z)), alpha = point_alpha, color = apal[[1]], position = position_jitter(), size = .5) +
geom_point(color = apal[[3]], alpha = 1) +
geom_line(color = apal[[3]]) +
geom_rug(data = unique(newdata, by = 'id'), aes(y = z_prime_id, x = NULL), color = apal[[4]], alpha = 1, sides = 'tr', length = unit(0.03, "npc")) +
geom_hline(yintercept = 0, color = apal[[1]]) +
coord_cartesian(ylim = c(-.025, .65), xlim = c(1, 10)) +
labs(y = 'FC correlation', x = 'Session') +
jftheme
}
max_points <- unlist(sea_schaefer400_7_withinnet[, .N, by = net1][, list(max = max(N))])
point_alphas <- sea_schaefer400_7_withinnet[, .N, by = net1][, p := 1 - N / max_points][, pomp := (p - min(p)) / (max(p) - min(p))][, alpha := .025 + pomp*.1]
font_add_google("Didact Gothic", "Didact Gothic")
showtext_auto()
for(i in 1:length(model_fits)){
model_fit <- model_fits[[i]]
network_name <- networks[[i]]
cat('\n\n### ', network_name, '{.tabset}\n\n')
cat('\n\n#### Plot\n\n')
point_alpha <- point_alphas[net1 == network_name, alpha]
print(plot_predictions(model_fit, point_alpha = point_alpha) + plot_percents(proportion_variance_tables[[i]]) +
plot_layout(ncol = 2, widths = c(4,1)))
cat('\n\n#### Table (%)\n\n')
print(knitr::kable(tidyr::spread(proportion_variance_tables[[i]], out_of, percent), digits = 1))
cat('\n\n#### Model Summary\n\n')
cat('\n```\n')
print(summary(model_fit))
cat('\n```\n')
}
| source | total | RE |
|---|---|---|
| ID | 0.9 | 14.1 |
| ID/Session | 5.7 | 85.9 |
| Residual | 93.4 | NA |
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 49086
Scaled residuals:
Min 1Q Median 3Q Max
-5.5418 -0.6779 -0.0802 0.5988 5.6291
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.004699 0.06855
id (Intercept) 0.000770 0.02775
Residual 0.076809 0.27714
Number of obs: 176897, groups: sess:id, 282; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.196427 0.006549 30
| source | total | RE |
|---|---|---|
| ID | 0.7 | 11.4 |
| ID/Session | 5.3 | 88.6 |
| Residual | 94.1 | NA |
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 182326.8
Scaled residuals:
Min 1Q Median 3Q Max
-4.7489 -0.6822 -0.0861 0.5939 6.3316
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.0045427 0.0674
id (Intercept) 0.0005854 0.0242
Residual 0.0813511 0.2852
Number of obs: 550247, groups: sess:id, 282; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.185095 0.005989 30.9
| source | total | RE |
|---|---|---|
| ID | 1.9 | 21.9 |
| ID/Session | 6.9 | 78.1 |
| Other RE | 0.0 | 0.0 |
| Residual | 91.2 | NA |
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 60958
Scaled residuals:
Min 1Q Median 3Q Max
-4.9158 -0.6941 -0.0934 0.5960 5.6157
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.006838 0.08269
id (Intercept) 0.001915 0.04376
Residual 0.090416 0.30069
Number of obs: 137824, groups: sess:id, 282; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.250244 0.009432 26.53
| source | total | RE |
|---|---|---|
| ID | 2.2 | 16.7 |
| ID/Session | 10.9 | 83.3 |
| Other RE | 0.0 | 0.0 |
| Residual | 86.9 | NA |
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: -1950.7
Scaled residuals:
Min 1Q Median 3Q Max
-4.0595 -0.6603 -0.0535 0.5948 5.3051
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.006863 0.08285
id (Intercept) 0.001377 0.03711
Residual 0.054545 0.23355
Number of obs: 39076, groups: sess:id, 267; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.189396 0.008577 22.08
| source | total | RE |
|---|---|---|
| ID | 1.3 | 14.7 |
| ID/Session | 7.5 | 85.3 |
| Other RE | 0.0 | 0.0 |
| Residual | 91.2 | NA |
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 27617.2
Scaled residuals:
Min 1Q Median 3Q Max
-5.9454 -0.6711 -0.0696 0.6003 6.6787
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.0057866 0.07607
id (Intercept) 0.0009995 0.03162
Residual 0.0702348 0.26502
Number of obs: 145733, groups: sess:id, 282; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.22049 0.00738 29.88
| source | total | RE |
|---|---|---|
| ID | 1.6 | 15.1 |
| ID/Session | 9.0 | 84.9 |
| Residual | 89.4 | NA |
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 132322.3
Scaled residuals:
Min 1Q Median 3Q Max
-6.2422 -0.6741 -0.1113 0.5556 7.7366
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.008211 0.09061
id (Intercept) 0.001455 0.03815
Residual 0.081652 0.28575
Number of obs: 393548, groups: sess:id, 282; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.277288 0.008834 31.39
| source | total | RE |
|---|---|---|
| ID | 0.6 | 13.7 |
| ID/Session | 4.0 | 86.3 |
| Residual | 95.3 | NA |
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 169633.5
Scaled residuals:
Min 1Q Median 3Q Max
-3.7405 -0.7070 -0.1085 0.6244 5.5554
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.0049931 0.07066
id (Intercept) 0.0007933 0.02817
Residual 0.1179362 0.34342
Number of obs: 240743, groups: sess:id, 282; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.267597 0.006692 39.99